function makenoise()
%
% function to make noise for stimulation. Noise can be filtered.
% Output file structure is for wavesamp in datac:
% file is ascii floating point, one entry per line
% sample rate is written on first line
% each block of npts (nominally 2048) entries are then read into
% the dac for successive trials
%
%fid = fopen('c:\pbm\datac\axo200b\noise02.txt', 'wt');
fsamp = 10000; % samples per second
fco = 100;		% cutoff frequency in Hz
wco = fco/(0.5*fsamp); % wco of 1 is for half of the sample rate, so set it like this...
%fprintf(fid, '%8.4f\n', 1000000/(2*fsamp)); % first entry is the sample rate. in usec/point (for 2 channels)
npts = 2^16;
xt=[0:1000/fsamp:(npts-1)*(1000/fsamp)]; % time base, in milliseconds
winl = 100;
w=hanning(winl*2);
w=w(1:length(w)/2); % just use half of it
yn(1)=0;
for j = 1:1 % make 100 independent samples of noise in this data set
   yn(2:npts)=randn(1,npts-1); % normally distributed random noise, 0 mean amplitude 1
   x=yn*0; %Allocate memory for x to speed up the next loop
%   for i=2:length(yn)
%      x(i)=(yn(i)+x(i-1))/3; % slightly correlated noise...
%   end
   %yn = x;
   % window the traces first
   yn(1:winl) = yn(1:winl).*w';
   n1 = length(yn) - winl + 1;
   n2 = length(yn);
   yn(n1:n2) = yn(n1:n2).*fliplr(w');
   %yn = x;
   if(wco < 1) % if wco is > 1 then this is not a filter!
%      [b, a] = ellip(4, 0.1, 60, wco)
[b, a] = butter(8, wco);
vsmo = filtfilt(b, a, yn); % filter all the traces..
   else
      vsmo = yn; % unfiltered version
   end;
   vsmo = vsmo - mean(vsmo);
   
   
   hf=findobj('Tag', 'plot');
   if(isempty(hf))
      hf = figure('Tag', 'plot');
   else
      clf(hf);
   end;
   [p, f] = pwelch(vsmo,  1024, fsamp);
   [px, fx] = pwelch(yn, 1024, fsamp);
   r = p./px; % this is the actual filtering that took place: the ratio of in to out...
   subplot('Position', [0.1 0.4 0.8 0.55]);
   %semilogx(f, r(:,1));
   loglog(f, p(:,1), 'g-', fx, px(:,1), 'k.');
   grid;
   subplot('Position', [0.1, 0.05, 0.8, 0.3]);
   plot(xt, vsmo);
%   fprintf(fid, '%8.3f\n', vsmo(1:2048));
end;
%fclose(fid);
return;
